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Abstract. A magic square is an n x n array of distinct positive integers whose sum along any 
row, column, or main diagonal is the same number. We compute the number of such squares for 
n — 4, as a function of either the magic sum or an upper bound on the entries. The previous record 
for both functions was the n — 3 case. Our methods are based on inside-out polytopes, i.e., the 
combination of hyperplane arrangements and Ehrhart's theory of lattice-point enumeration. 



1. It's A Kind Of Magic 

A magic square is an n x n array of distinct positive integers whose sum along any row, column, or 
main diagonal is the same number, the magic sum. The history of magic squares is well documented, 
see, e.g., [8, 9, 21]. The contents of a magic square have varied with time and writer; usually 
they have been the first n 2 consecutive positive integers, but often any arithmetic sequence and 
sometimes fairly arbitrary numbers. The fixed ideas are that they are integers, positive, and 
distinct. 

In the last century mathematicians took an interest in results about the number of squares with 
a fixed magic sum, but with simplifications: diagonal sums were often omitted and the fundamental 
requirement of distinctness was almost invariably neglected [1, 3, 11, 16]. For example, classical 
formulas of MacMahon [15] include 

t + 3\ ft + 2 

the number of 3 x 3 squares with (not necessarily distinct) nonnegative integer entries that sum to 
t along any row and column, and 

-t 2 + -t+l if 3|t, 

the number of such squares in which the two main diagonals also sum to t. The papers [6, 20] form, 
to the best of our knowledge, the beginning of a theory that tackles counting problems related to 
magic squares with the distinctness of the entries enforced. 

Our goal is to show that the ideas in [6] can be used to compute the number of 4 x 4 magic squares 
(with distinct entries), as a function of a parameter that is either the magic sum or an upper bound 
on the entries of the square. To be precise, we define the affine magic counting function a n (t) to be 
the number of all n x n matrices consisting of distinct positive integers whose sum along any row, 
column, or main diagonal is the same number t. The cubical magic counting function c n (t) is the 
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number ofnxn matrices whose entries are distinct positive integers less than t, whose sums along 
any row, column, or main diagonal are equal. The previous record consisted of the 3x3 counting 
functions [6, 20] (see also [5] for the computational implementation of [6]) 
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Both of these functions are quasipolynomials, i.e., of the form c<i(t) t d + Cd-i(t) t^ 1 + • • • + co(t), 
where Co, c±, . . . , are periodic functions. It follows from [6] that a n {t) and c n {t) are always 
quasipolynomials in t; we will outline the basic arguments in Section 2. 

A handy, compact way of representing a quasipolynomial q(t) is through its generating function 
Q( z ) := St>o9(*) z< - I* i s n °t nar d to prove (see, e.g., [4, 17]) that Q(z) is a rational function 
with poles at pth roots of unity, where p is a common period of the coefficient functions of q{t). As 
examples we give the rational generating functions for 03 (t) and cs(t): 

8x 15 (2x 3 + l) 



(l-x 3 )(l-x 6 )(l 



and 
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Our computational approach (which differs from both [6] and [20]) verified the results for 03 (t) 
and cs(t). Both papers [6, 20] commented that the 4x4 case seems computationally infeasi- 
ble, and we present some reasons for this assessment in Section 3; one of the reasons is that 
the rational generating functions for 04 (t) and C4(i) take several pages to write down (in re- 
duced form). Nevertheless we were able to compute A^(z) and C/±{z); since we did not want 
to waste paper, the results are posted online at math.sfsu.edu/beck/papers/affmagic4.html 
and math. sfsu. edu/beck/papers/cubmagic4. html. 

2. Enter Geometry 

2 

We start with the affine case. One treats a square with magic sum t as an integer vector x £ Z n 
confined to the affine subspace ts, where 



x G 



all line sums 



equal 1 j . 
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From a geometric point, a square with magic sum t is an integer point in the open dilated polytope 
tV = M>o (Its, except that we need to require that the entries are distinct; so our integer point has 
to lie outside of the hyperplane arrangement 

H = {xij = x k i : 1 < i,j,k,l < n, / (k,l)} . 

The pair (V, 7i) is referred to as an inside-out polytope [7]. One can now use Mobius inversion on this 
hyperplane arrangement (the computational approach in [6] ; one appealing feature of this approach 
is that one computes generating functions for polytopes of progressively smaller dimension) or 
view (M>o fl ts) \ (J H as a disjoint union of open polytopes (our approach, which we call region 
enumeration). Either way, Ehrhart's fundamental results on integer-point enumeration in polytopes 
[4, 10] imply that a n (t) is a sum/difference of Ehrhart quasipolynomials — lattice-point enumerators 
of rational polytopes. 

The cubical case is very similar, except that we replace the subspace by 

s = jx G M™ 2 : all line sums equal j 

and the polytope by (0, 1)™ Pi s. The hyperplane arrangement, which enforces distinct entries in 
the square, is the same as in the affine case. 

3. Computations 

For the affine case (V = R>o fl s), the present computational approach treats tV \ \JH as a 
disjoint union of open polytopes (the regions of the inside-out polytope (V,7i)) and computes 
A^{z) by summing the Ehrhart generating functions of all of the polytopes. The computation of 
C&{z) in the cubical case is similar. The challenge in this approach is computing the regions and 
summing the resulting rational functions efficiently. Our algorithm for computing the regions of 
an inside-out polytope is fully described in [13] and runs in polynomial time for polytopes of fixed 
dimension. 

To briefly outline the process, we observe that the removal of a single transverse hyperplane 
(i.e., one that intersects the relative interior of V) from an open polytope splits the polytope into 
two new open polytopes. Removing a non-transverse hyperplane simply results in the same open 
polytope. To enumerate the regions of an inside-out polytope, we begin by splitting the original 
polytope with the first hyperplane. We then proceed recursively, splitting the newly generated 
polytope(s) with the next hyperplane, and so on, until all the hyperplanes in the arrangement are 
exhausted. To determine whether a hyperplane is transverse to a given polytope, we engage a linear 
solver, using a normal vector to the hyperplane as the objective function, and the polytope as the 
feasible region. Computation of each region's Ehrhart generating function is accomplished using 
Barvinok's algorithm [2]. 

Aside from the running time of Barvinok's algorithm (which is polynomial in fixed dimension), 
the computational complexity of our approach has a similar order of magnitude to Mobius inversion, 
namely Yli=i (") where d = dimV and n = \H\. Implemented as a depth-first traversal of a binary 
tree [13], memory requirements are linearly proportional to \H\. Despite the disadvantage of needing 
to compute Ehrhart generating functions for maximal-dimensional polytopes (as opposed to the 
Mobius inversion approach), region enumeration lends itself well to parallel computation. At a 
predetermined level of recursion, we send a description of each open polytope, together with the 
remaining hyperplanes, to separate processors and accumulate the results. 

We have implemented the region enumeration algorithm using C++, combining the polytope 
libraries barvinok [19], polylib [14], and cdd [12] with the computer algebra system GP/Pari [18] 
to form a suite of applications and a small library capable of operations on inside-out polytopes, 
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collectively referred to as IDP [13] (http://iop.sourceforge.net/). As far as we are aware, 
IOP is the first software implementation specifically capable of computing generating functions 
of inside-out polytopes. The library includes the ability to convert an inside-out polytope to an 
isomorphic, lattice-point equivalent problem with full dimension in the ambient space, a feature 
that significantly reduces the time required to generate the regions. 

To compute the number of 4 x 4 magic squares, we have reduced the scope of each problem 
by taking advantage of the 32 symmetries described in [13]. The reduced problems were fed to 
a 40-node Dell PowerEdge cluster. Computing the generating functions A±{t) and C±(t) each 
took less than a day. Each reduced problem required the enumeration of 3211412 polytopes and 
the computation, summation, and simplification of an equal number of rational generating func- 
tions. In lowest terms, the rational generating functions A±{i) and C±{t) respectively occupy ap- 
proximately 280k and 25k of disk space with numerator /denominator polynomials of degree 2900 
(affine) and 533 (cubical) and coefficients exceeding 10 42 . Files containing the results are available 
at math. sfsu.edu/beck/papers/ [affmagic4.html , cubmagic4 .html] . The extreme nature of the 
simplified results justifies the skepticism voiced in [6] and [20] about the feasibility of computing 
counting functions for 4x4 magic squares. 

During development of IOP we observed phenomena that may be worthy of future investigation. 
The degree to which rational-function arithmetic affected computation time was surprising. In 
a test run for the affine case, a 1.4 gHz home computer successfully enumerated the regions for 
a portion of the reduced problem [13, Equation 6.16] and output individual simplified Ehrhart 
generating functions to a file in approximately a week's continuous run-time. However, the attempt 
to generate a simplified sum proved intractable, with a projected run-time of well over three months. 
This is considerably more gHz-hours than were actually used by the Dell cluster to tackle the entire 
reduced problem, suggesting that associative grouping of the rational generating functions may 
play a significant role in the run-time. In particular, we suspect that grouping rational functions in 
expressions corresponding to an entire inside-out polytope may result in smaller simplified forms 
(hence, faster computation time) than expressions grouped arbitrarily. 
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